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ABSTRACT 


The  Pseudo  Wigner-Ville  Distribution  is  a  signal  transformation  of  an  input  time  signal 
into  a  joint  time-frequency  domain,  that  provides  an  excellent  characterization  of  an  input 
signal  as  well  as  it’s  respective  energy  content.  A  smoothing  window  and  energy  sensitivity 
analysis  of  the  discretized  Pseudo  Wigner-Ville  Distribution  is  presented  along  with  a 
machinery  monitoring  application  using  the  resulting  Wigner-Ville  energy  of  the  sampled 
signal.  The  ability  to  apply  the  Pseudo  Wigner-Ville  Distribution  to  both  steady-state  and 
transient  machinery  may  enable  the  monitoring  and  diagnostics  of  virtually  any  piece  of 
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I.  INTRODUCTION 


The  recent  world  events  have  caused  drastic  changes  in  many  areas  and  made  the  future 
that  much  more  unpredictable.  With  this  high  degree  of  uncertainty  and  drastic  financial 
cutbacks  as  evidenced  in  the  Defense  budget,  the  emphasis  of  technological  advance  and 
dominance  is  all  the  more  heightened.  With  the  requirement  to  ”do  more  with  less",  it  is 
imperative  to  find  ways  to  save  money.  It  is  with  this  mind-set  that  better  ways  to  prolong 
the  life  of  operating  machinery  and  defray  unneeded  maintenance  costs  be  investigated.  For 
many  years,  most  transient  machinery  has  gone  unmonitored  and  removed  at  periodic  time 
intervals  for  refurbishment  regardless  of  the  current  machinery  condition.  Developing  a 
method  to  monitor  and  diagnose  these  machines  of  a  transient  nature  will  provide  a  way  to 
determine  when  overhaul  or  replacement  is  needed. 

Vibrational  data  has  long  been  used  as  the  tool  to  measure  a  machine’s  health.  Early 
methods  of  monitoring  machinery  had  been  for  the  experienced  personnel  to  listen  or  touch 
the  machine  to  detect  if  a  fault  exists,  however,  this  method  is  dependent  on  human  senses 
and  the  same  individual  conducting  the  monitoring.  It  was  for  this  reason  that  vibration  levels 
as  monitored  by  a  dynamic  signal-analyzer  are  now  used  for  machinery  health  monitoring. 
The  vibration  levels  of  a  good  condition  establish  a  baseline  to  allow  a  means  of  comparison 
when  subsequent  measurements  are  taken.  These  dynamic  signal  analyzers  are  time 
independent  and  do  not  show  the  time  dependency  required  by  transient  machinery.  It  is 
because  of  this  limitation,  that  a  three-dimensional  time,  f requency  and  energy  representation 
such  as  the  Pseudo  Wigner-Ville  Distribution  is  being  investigated  as  a  method  to  monitor 
transient  machinery. 
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A.  PREVIOUS  WORK  AND  ACCOMPLISHMENTS 

This  thesis  is  a  follow-on  work  which  has  been  preceded  by  LT  G.  Rossano  and  LCDR 
J.  Calahorrano,  both  in  conjunction  with  Prof.  Y.S.  Shin.  Rossano  established  the  initial  form 
of  the  computer  program  that  calculates  the  discretized  Pseudo  Wigner-Ville  Distribution. 
He  evaluated  the  effect  that  different  signals  and  procedures  had  on  the  Distribution  along 
with  characterizing  an  actual  pump  signal.  Calahorrano  made  some  further  improvements  to 
the  computer  program  and  then  applied  this  improved  version  of  the  computer  program  to 
characterize  some  classified  signals.  Both  of  these  works  have  been  a  great  benefit  in  dealing 
with  characterizing  signals,  but  have  not  investigated  the  smoothing  method  and  energy  or 
amplitude  of  the  discretized  Pseudo  Wigner-Ville  Distribution. 

B.  OBJECTIVES  OF  CURRENT  RESEARCH 

The  Pseudo  Wigner-Ville  Distribution  provides  an  outstanding  time-frequency  domain 
spectra  representation  of  an  input  signal  as  shown  through  the  works  among  others  of  Rossano 
and  Calahorrano  here  at  the  Naval  Postgraduate  School.  The  purpose  of  this  research  is: 

•  To  review  and  modify  the  computer  program  developed  here  at  the  Naval  Postgraduate 
School  to  support  the  smoothing  window  and  energy  analysis  work. 

•  To  conduct  a  smoothing  window  sensitivity  analysis. 

•  To  investigate  the  energy  content  of  the  discretized  Pseudo  Wigner-Ville  Distribution. 

•  To  apply  the  developed  smoothing  window  and  energy  principles  to  an  actual 
machinery  monitoring  application. 

•  To  investigate  the  feasibility  of  using  the  Pseudo  Wigner-Ville  Distribution  energy  as 
a  pre-processed  input  to  a  Neural  Network. 
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11.  THE  PSEUDO  WIGNER-VILLE  DISTRIBUTION 


A.  THEORETICAL  BACKGROUND 
1.  Evolution 

The  current  day  Pseudo  Wigner-Ville  Distribution  is  a  three  dimensional  (time, 
frequency,  and  energy)  representation  of  a  signal  that  is  particularly  well  suited  for  analysis 
of  non-staticnary  signals.  Eugene  Wigner  [Ref.  1]  first  introduced  the  Wigner  distribution  in 
1932  to  study  the  problem  of  statistical  equilibrium  in  the  area  of  quantum  mechanics.  His 
work  was  furthered  in  1948  by  J.  Ville  [Ref.  2],  who  used  the  Wigner  distribution  in  the  area 
of  signal  analysis.  A  major  contribution  of  Ville’s,  is  the  use  of  analytic  signals  as  an  input 
to  the  distribution  vice  the  customary  real  signal.  The  advantages  of  using  an  analytic  signal 
in  the  distribution  are  two-fold  [Ref.  3].  First,  the  distribution  of  a  real  signal  results  in  a 
symmetrical  spectrum  with  only  half  of  the  result  containing  useful  information,  while  the 
use  of  an  analytic  signal  avoids  this  negative  frequency  redundancy.  Secondly,  by  accounting 
for  only  the  positive  frequencies  it  satisfies  both  the  practical  and  the  mathematical 
completeness  of  the  problem.  The  third  part  of  the  Pseudo  Wigner-Ville  Distribution  (ie. 
Pseudo)  ,  arises  from  the  need  to  apply  a  smoothing  window  to  the  resulting  distribution. 
This  is  a  result  of  the  presence  of  cross  terms  that  arise  from  multiple  component  signals  and 
will  be  covered  in  depth  later. 

The  works  of  T.A.C.M.  Claasen  and  W.F.G.  Mecklenbrauker  [Refs.  4,5,  and  6]  in 
1980  have  paved  the  way  for  many  recent  applications  of  the  Pseudo  Wigner-Ville 
Distribution.  Their  three  part  paper  is  an  all  encompassing  work  that  has  served  to  highlight 
the  capabilities  of  the  Pseudo  Wigner-Ville  Distribution  and  allow  for  greater  knowledge  of 
the  subject  with  an  emphasis  on  signal  processing  and  analysis. 
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A  considerable  amount  of  recent  work  using  the  Pseudo  Wigner-Ville  Distribution 
have  come  in  the  areas  of  optics  [Refs.  7,8,  and  9]  and  speech  [Ref.  10  and  11].  Additionally, 
T.J.  Wahl  and  J.S.  Bolton  of  Purdue  University  [Ref.  12]  have  used  the  Pseudo  Wigner-Ville 
Distribution  to  analyze  structural  impulse  responses.  Two  recent  papers  have  investigated  the 
use  of  the  Wigner-Ville  Distribution  in  the  area  of  machinery  monitoring.  Flandrin  et.  al. 
[Ref.  13]  proposed  the  Wigner-Ville  as  a  means  to  confirm  a  machinery  diagnosis  in  the 
specific  application  of  a  Pressurized  Water  Reactor  power  plant  and  lastly,  Forrester  [Ref.  14] 
has  investigated  the  Wigner-Ville  Distribution  as  a  method  for  fault  detection  in  gears.  As 
evidenced  by  the  cited  examples,  the  capability  and  adaptability  of  the  Wigner-Ville 
Distribuiion  is  enormous  with  applications  of  it's  use  rapidly  expanding. 

2.  Function  Definition 

The  Wigner-Ville  Distribution  is  a  transformation  of  a  signal  into  the  time- 
frequency  domain.  By  definition,  the  cross  Wigner-Ville  Distribution  is  defined  as; 

dx 

where;  r(t)  =  a  complex  time  history 
s(t)  =  a  complex  time  history 
t  =  time 
w  =  frequency 
•  =  complex  conjugate 

Similarly,  the  auto  Wigner-Ville  Distribution  is  defined  as; 

m 

VfDF^^~  re-J“^  r(t+-^)  dx 
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Since  the  concentration  of  this  research  is  with  the  representation  of  a  single  input  signal,  the 
auto  Wigner-Ville  Distribution  will  be  used.  This  expression  is  essentially  a  Fourier  transform 
of  the  auto-correlation  of  a  signal,  which  may  be  thought  of  as  a  three  dimensional  spectral 
density  function. 

3.  Distribution  Properties 

Thwe^^are  several  properties  of  the  Distribution  that  are  worth  noting.  The 
properties  listed  below  have  been  shown  in  many  works  concerning  the  Wigner-Ville 
Distribution.  These  properties  will  be  shown  in  the  form  of  the  auto  Wigner-Ville 
Distribution  [Ref.  15]. 

Time  and  frequency  shift  properties  follow: 

•  A  time  shift  in  the  signal  corresponds  to  a  time  shift  of  the  Distribution: 

•  A  frequency  shift  in  the  signal  corresponds  to  a  frequency  shift  in  the  Distribution: 

•  Both  a  frequency  and  time  shift  in  the  signal  corresponds  to  the  same  in  the 
Distribution: 


Since  our  concern  is  in  the  area  of  machinery  monitoring  and  diagnostics,  the  aforementioned 
properties  are  a  nessesity  in  the  development  of  a  suitable  monitoring  program. 
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The  next  three  properties  provide  information  concerning  the  energy  contained 
within  the  Wigner-Ville  Distribution. 

•  The  integration  of  the  Distribution  over  the  frequency  domain  provides  the 
instantaneous  signal  power: 


-  I  s(t)  P 


•  The  integration  of  the  Distribution  over  the  time  domain  provides  the  energy  density 
spectrum: 


f  -  1  5(w)  P 


•  The  integration  over  both  the  time  and  frequency  domain  provides  the  total  energy  of 
the  input  signal: 


-~ff  WDF^^g(t,Ci>)  dtdo>  -  I  I  s  I  P 

Again,  as  with  the  time  and  frequency  properties,  the  energy  contained  in  the  Wigner-Ville 
Distribution  must  be  definable  in  order  to  use  this  in  a  machinery  monitoring  program. 
Additional  properties  may  be  found  in  Reference  4. 

4.  Cross  Terms  Within  The  Distribution 

A  particularly  troublesome  characteristic  of  the  Wigner-Ville  Distribution  is  the 
presence  of  cross  terms.  These  cross  terms  result  in  the  Distribution  taking  on  negative  values 
as  shown  later  in  Figure  2.  Since  the  properties  above  have  shown  that  the  Wigner-Ville 
Distribution  represents  an  energy  value,  these  negative  values  raise  an  area  of  concern.  Jones 
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and  Parks  [Ref.  16]  show  that  these  cross  terms  are  essentially  the  interference  of  multiple 
component  signals  that  occur  when  there  is  a  non-zero  cross-correlation  within  the  signal. 
The  equation  listed  below  shows  the  cross-term  that  results  from  a  signal  that  consists  of  two 
components. 

WDF^,y(t.U>)  -  +  2RE[WDF^y(t,Ui)]  +  WDFy(t,Ui) 


As  shown,  the  cross-term  is  represented  by  the  center  term  while  the  other  terns  represent 
their  respective  signal  components.  There  are  several  methods  for  eliminating  these  terms  and 
providing  a  more  presentable  Distribution  and  there  has  been  a  considerable  amount  of 
research  in  the  are  of  an  optimum  smoothing  method  for  the  Wigner-Ville  Distribution  as 
indicated  by  References  3,4,17,  and  18.  The  method  chosen  for  smoothing  the  Wigner-Ville 
Distribution  in  this  work  is  a  post  processing,  sliding  Gaussian  window  function.  What  is 
meant  by  post  processing,  is  that  the  distribution  is  smoothed  after  calculation  of  the  Wigner- 
Ville  Distribution. 

5.  Gaussian  Smoothing  Function 

The  Gaussian  window  function  was  first  proposed  mathematically  as  a  means  to 
smooth  the  Wigner-Ville  Distribution  by  N.D.  Cartwright  [Ref.  19].  The  continuous  form  of 
the  Gaussian  window  function  follows; 


G{  t,  u) 


2710 


Wahl  and  Bolton  in  Reference  12  have  furtherd  this  work  and  developed  the  sampled  form 
of  the  window  function  that  is  used  in  our  discretized  form  of  the  Pseudo  Wigner-Ville 
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Distribution  for  smoothing.  The  effect  of  various  size  smoothing  windows  will  be  presented 
in  greater  depth  in  Chapter  III. 

B.  WIGNER-VILLE  DISTRIBUTION  COMPUTER  PROGRAM 

The  computer  program  listed  in  the  Appendix  is  a  method  that  was  developed  here  at 
the  Naval  Postgraduate  School  for  calculating  the  discretized  Pseudo  Wigner-Ville 
Distribution.  As  mentioned  earlier,  the  program  that  computes  the  <!iscretized  Pseudo 
Wigner-Ville  Distribution  is  largely  the  work  of  G.  Rossano  [Ref.  15],  who  initially 
investigated  the  use  of  the  Wigner-Ville  Distribution  as  a  tool  for  machinery  monitoring.  For 
the  purposes  of  this  work,  the  program  has  been  thoroughly  reviewed  and  rewritten  to  allow 
for  ease  of  readibility  and  modification.  The  current  version  enables  the  user  to  look  at  not 
only  the  resulting  smoothed  version  of  the  Pseudo  Distribution,  but  also  the  unsmoothed 
Wigner-Ville  Distribution.  The  second  major  modification  was  to  program  the  smoothing 
window  subroutine  to  allow  for  the  use  of  different  size  windows  in  smoothing  the 
Distribution.  These  changes  have  proven  invaluable  for  evaluating  the  effect  of  various  sized 
smoothing  windows  and  in  conducting  the  energy  sensitivity  analysis.  Presented  below  are 
the  two  major  equations  that  comprise  the  calculation  of  the  discretized  Pseudo  Wigner-Ville 
Distribution  as  performed  in  the  program  listed  in  the  Appendix. 

1.  Discrete  Time  Wigner-Ville  Distribution 

As  mentioned  earlier,  the  Wigner-Ville  Distribution  is  essentially  a  Fourier 
transform  of  an  auto-correlation  function.  There  have  been  a  couple  of  equations  developed 
to  represent  the  discrete  Wigner-Ville  Distribution.  The  version  that  has  been  used  in  the 
development  of  the  program  listed  in  the  Appendix,  is  one  that  was  developed  in  Reference 
4  by  Classen  and  Mecklenbrauker  and  and  later  used  in  the  work  of  Wahl  and  Bolton  [Ref. 
12]  and  Rossano  [Ref.  IS].  The  only  difference  of  note  between  the  various  discrete  versions 
is  essentially  a  difference  in  constants.  All  of  the  versions  are  basically  a  Fourier  transform 
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of  an  auto-correlation  function.  Shown  below  is  the  discretized  equation  that  has  been  used 
for  programming  purposes  in  this  work: 

N 

^(iAto, jAt)  -  y^Re{FFT[2At  CORR(i)]) 

pi 


where:  Re  »  Real  part 

Corr(i)  >  complex  auto-correlation  of  the  signal  s(t) 
FFT  *  Fast  Fourier  Transform 
At  *  sample  time 


2.  Discretized  Gaussian  Window  Function 

Due  to  the  need  to  deemphasis  the  resulting  cross-terms  within  the  Wigner-Ville 
Distribution,  a  suitable  method  for  smoothing  has  been  obtained.  The  sampled  equation  of 
the  selected  smoothing  window  is  shown  below  [Ref.  12]: 


where: 


-2M<m<2M  & 

=  MAt  Si 


1 
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-2N<n^N 
=  NAw 


-( 


e 


(aiAt)* 

2o\ 


(nAw)  ^  , 
2oi 
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III.  SMOOTHING  WINDOW  ANALYSIS 


As  shown  above,  the  smoothing  function  or  window  that  is  used  in  this  work  is  a 
Gaussian  window  that  is  applied  to  the  unsmoothed  Wigner-Ville  Distribution.  Figure  1  is  the 
raw  time  input  signal  of  a  100  &  400  Hz  sine  wave  that  will  be  processed  through  the  Wigner- 
Ville  Distribution.  Figure  2  is  the  unsmoothed  Distribution  that  results  from  this  combined 
100  and  400  Hz  sine  wave.  As  can  easily  be  seen  in  Figure  2,  the  resulting  cross  terms  are 
very  dominant  and  need  to  be  deemphasized.  This  Gaussian  smoothing  function  is  applied 
after  the  Wigner-Ville  Distribution  has  been  calculated  (ie.  post  processing),  which  allows  for 
visual  inspection  of  both  the  unsmoothed  Wigner-Ville  Distribution  and  the  effect  that 
different  sized  windows  may  produce  on  the  smoothed  Distribution  (ie.  Pseudo).  For 
example.  Figure  3  is  the  smoothed  result  of  Figure  2  after  applying  a  13  by  13  sized  window. 
The  next  few  sections  will  discuss  some  of  the  characteristics  of  this  post-processing 
smoothing  window  and  what  may  be  gained  or  lost  through  varying  the  size  of  the  applied 
window. 
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100  0  400  HZ  SINE  WAVES 
8IQNAL  LENGTH  =  0.250  SBC  /  AMPS  =  1.0 


(8M0A)  ganindiAiv 


Figure  1.  100  &  400  Hz  Sine  Wave  Input 
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Figure  3.  Pseudo  Wigner-Ville  Distribution  (13  X  13  Window) 


A.  SIZE  CRITERION 


As  shown  by  N.D.  Cartwright  [Ref.  19]  for  the  continuous  form  of  the  Wigner-Ville 
Distribution,  the  size  of  the  smoothing  window  may  be  selected  such  that  a  positive 
Distribution  will  result.  Following,  are  these  criterion  that  have  been  mathematically  proven 
for  the  continuous  case  to  assure  a  positive  distribution; 


where: 

a  j  -  AfA  t  o^,  -  iVA(i> 

By  chosing  M  and  N  such  that  this  criterion  is  met,  the  resulting  smoothed  discretized 
Wigner-Ville  Distribution  does  not  necessarily  result  in  all  of  the  values  being  positive.  This 
may  be  seen  in  Figure  4,  where  the  window  that  was  used  was  again  a  13  by  13  (ie.  N  =  13 
&  M  *  13)  sized  window.  However,  it  does  deemphasize  the  large  cross-terms  to  the  point 
that  they  are  nearly  zero  and  yields  a  very  presentable  Pseudo  Wigner-Ville  Distribution.  The 
ability  to  deemphasize  these  cross  terms  is  a  necessity  if  the  desire  is  to  evalute  the  energy 
change  within  a  signal.  For  the  remainder  of  this  work,  the  window  size  selected  (ie.  N  &  M, 
N  X  M,  or  N  by  M)  will  be  such  that  this  criterion  is  met  as  closely  as  possible. 
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Smooth^ 


B.  TIME  AND  FREQUENCY  DEPENDANCE 

Since  the  criterion  shown  previously  are  time  and  frequency  dependent,  there  must  be 
a  relation  oetween  the  two  that  can  be  reduced  to  one  parameter.  The  customary  time  and 
frequency  relations  of  signal  analysis  are  slightly  changed  in  the  context  of  the  Wigner-Ville 
Distribution.  Shown  below  are  the  time  and  frequency  relations  when  working  with  the 
Wigner-Ville  Distribution. 

;  Af  -  ^  -  2xDPxAf 

where:  tn^^x  “  signal  length  &  f„j^  =  maximum  frequency 

At  =  sampling  time  &  Af  =  frequency  resolution 
DP  =  number  of  data  points 


The  definitions  of  the  size  criterion  discussed  earlier  for  the  smoothing  window  follow; 

2  2  1 
o\ai  ^ 

4 

and  since; 

a\  -  oi  -  (A7Aa>)2 

and  the  time  and  frequency  are  related  in  the  Wigner-Ville  Distribution  through  the  number 
of  data  points  as  shown; 


Aci) 


7C 

2  DP  At 
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Making  this  substitution  into  the  size  criterion,  yields  a  very  simple  result  for  determining  the 
size  of  the  smoothing  window  required. 

MxN  i.  — 
n 

In  looking  at  this  result  it  seems  as  if  the  window  is  completely  independent  of  time  and 
frequency,  but  it  must  be  rememberd  that  M  and  N  are  parameters  that  specify  the  lengths 
of  the  smoothing  window  in  the  time  and  frequency  directions  respectively. 

C.  RELATION  TO  NUMBER  OF  DATA  POINTS 

Another  area  of  interest  is  the  relation  of  the  window  size  to  the  number  of  data  points 
in  the  sampled  signal.  The  best  way  to  show  this  relation  is  through  an  example.  For 
instance,  a  signal  of  a  length  4.0  seconds  that  is  sampled  using  512  data  points,  will  result  in 
a  At  *  7.8125  msec  and  Af  =  0.0625  Hz  as  calculated  as*ng  the  constraints  of  time  and 
frequency  as  defined  by  the  Wigner-Ville  i  -usiribution.  Since  the  example  has  used  512  data 
points,  the  previous  section  has  shovn  that  chasing  N  =  10  and  M  =  16  will  approximately 
satisf y  the  size  criterion  limitations.  Additionaly,  by  definition  of  the  smoothing  window,  the 
area  of  coverage  will  be: 

-2MiW^2M  -2mn^2N 

This  means  that  for  the  example  shown,  the  coverage  of  the  window  will  be  40  steps  in  the 
frequency  domain  and  64  steps  in  the  time  domain. 
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40xAjf-  2.5  64xAt  -  0.5  AREA  25 

where:  Af  =  0.0625  Hz  &  At  «  7.8125  msec 

The  AREA  »  1.25  is  a  unitless  quantity  that  defines  the  area  that  the  smoothing  window 
encompasses  while  smoothing  one  point  in  the  Distribution.  The  unitless  dimension  is  the 
result  of  the  time  being  in  seconds  and  the  frequency  in  hertz.  For  comparison  purposes,  the 
total  area  of  the  distribution  is  the  product  of  f„,.,  and  t„.„  which  is  equal  to  256.0,  resulting 
in  the  window  smoothing  over  0.488%  of  the  total  Distribution,  while  smoothing  one  point. 

Now  consider  another  example  that  uses  1024  data  points.  As  shown  earlier,  the 
product  of  M  and  N  must  be  greater  than  or  equal  to  the  number  of  data  points  divided  by 
pi  in  order  to  meet  the  window  size  criterion  as  established  by  Reference  12.  With  respect  to 
the  earlier  example,  the  size  of  the  window  (ie.  N  times  M)  must  increase  two  times. 
Therefore,  a  selection  of  N  =  10  and  M  =  32  will  be  used  for  this  example.  These  values  are 
very  close  to  the  required  criterion  and  will  suffice  for  this  illustrative  example.  The  signal 
length  (t„„)  will  be  2.0  seconds  resulting  in  f^^  ■  256.25  Hz,  At  =  1.953  msec,  and  Af  = 
0.125  Hz.  In  relation  to  the  area  of  coverage  of  the  window,  the  following  results  apply: 

40xAf-5.0  128xAt  -  0.249984  AREA  -1.25 

where:  Af  =  0.125  Hz  &  At  =  1.953  msec 

This  is  the  same  result  as  the  512  data  point  case,  with  the  following  exception.  The  total  area 
of  the  Distribution  (t„„  x  f^^  =  512.0)  has  doubled,  while  the  area  of  the  smoothing  window 
coverage  has  remained  the  same.  This  results  in  the  window  only  covering  0.244%  of  the 
entire  Distribution  when  smoothing  a  single  point  within  the  Distribution.  In  conclusion,  the 
greater  the  number  of  sampled  data  points  within  the  processed  raw  time  signal,  the  smaller 
the  area  the  smoothing  window  will  encompass  while  smoothing  the  Distribution.  An 
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important  trade-off  to  consider,  however,  is  the  increase  in  computer  time  required  to  run 
the  additional  number  of  data  points. 

D.  TIME/FREQUENCY  RESOLUTION 

Since  the  smoothing  window  used  in  this  work  gives  the  ability  to  vary  the  time  and 
frequency  size  of  the  window,  the  following  examples  show  the  results  that  can  be  obtained. 
Figure  S,  is  the  input  signal  that  is  sampled  using  SI 2  data  points.  It  is  a  400  Hz  sine  wave 
that  has  been  computer  generated  with  an  amplitude  of  2.0.  In  order  to  be  consistent  with  the 
previous  work,  the  size  criterion  established  earlier  will  be  used  as  a  guide  for  the  selection 
of  M  and  N  in  the  smoothing  window.  Figure  6,  is  the  result  of  a  window  in  which  M  and 
N  are  both  equal  to  13.  Since  this  smoothing  window  does  not  require  M  and  N  to  be  equal, 
varying  these  parameters  yielded  some  intersting  results.  Figure  7,  shows  the  result  of  the 
input  signal  being  smoothed  with  a  35  by  S  window.  This  means  that  N  =  35  and  M  =  5,  or 
in  other  words,  the  window  is  longer  in  the  frequency  domain  than  it  is  in  the  time  domain. 
The  resulting  distribution  has  an  amplitude  that  is  much  more  constant  throughout  the  entire 
time  length,  but  the  amplitude  itself  has  been  reduced  considerably  as  compared  to  the  13  by 
13  example.  Figure  8  is  just  a  different  view  of  this  example,  showing  the  loss  of  frequency 
resolution  with  a  35  by  5  window.  Just  the  opposite  case  to  this  example  is  where  the  window 
used  for  smoothing  is  a  5  by  35  window.  In  this  case,  the  window  is  longer  in  the  time 
domain  than  the  frequency  domain  and  just  the  opposite  results  are  obtained.  Figure  9,  shows 
the  result  of  this  window,  while  Figure  10  gives  a  view  showing  the  degree  of  frequency 
resolution  that  can  be  obtained.  This  can  be  very  beneficial  in  the  monitoring  of  machinery, 
since  in  many  cases  the  ability  to  differentiate  between  signals  that  are  very  close  to  one 
another  is  paramount. 
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400  HZ  SINE  WAVE 

SIONAL  LENGTH  s  0.260  SEC  /  AMP  a  2.0 
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Figure  6.  Pseudo  Wigner-Ville  Distribution  (13  X  13  Window) 
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Figure  7.  Pseudo  Wigner-Ville  Distribution  (35  X  5  Window) 


Amp.  (V“2/{H2*sec)) 


400  Hi  SINE  WAVE 
FREQUENCY  RESOLUTION  ;  35  X  S 
AMPUTUDE  =  2.0 


Frequency  (Hz) 


Figure  8.  Frequency  Resolution  (3S  X  S  Window) 
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Figure  9.  Pseudo  Wigner-Ville  Distribution  (5  X  35  Window) 
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Amp.  (V“2/(Hz*sec)) 


400  Hz  SINE  WAVE 
FREQUENCY  RESOLUTION  ;  5  X  3S 
AMPLITUDE  =  2.0 


250.0  500.0  750.0 

Frequency  (Hz) 


Figure  10.  Frequency  Resolution  (5  X  35  Window) 


E.  TRANSITIONAL  PEAK  PHENOMENON 

In  works  that  have  given  examples  of  the  Wigner-Ville  Distribution  for  transient 
signals,  the  point  at  which  the  signal’s  frequency  is  changing  from  increasing  to  decreasing 
or  vice  versa  has  a  large  amplitude  increase  in  the  Wigner-Ville  Distribution.  For  example, 
Rossano  [Ref.  IS]  shows  an  artificially  generated  linear  chirp  and  also  an  actual  pump 
signature  that  has  what  appears  to  be  a  large  concentration  of  energy  at  the  transition  point. 
In  order  to  use  the  Wigner-Ville  Distribution  for  transient  machinery  monitoring,  there  must 
be  an  explanation  or  way  to  quantif y  this  "ghost  peak"  that  is  occuring  at  the  transition  point. 
To  date,  there  has  not  been  an  explanation  to  this  phenomenon,  however,  by  looking  at  the 
unsmoothed  Wigner-Ville  Distribution,  the  reason  for  this  is  easily  seen.  Figure  1 1  is  a  linear 
chirp  sine  wave  that  is  used  as  the  input  to  the  Wigner-Ville  Distribution.  Figure  12,  shows 
the  resulting  unsmoothed  Wigner-Ville  Distribution  of  this  linear  chirp  and  as  discussed  and 
shown  earlier,  the  cross-terms  can  be  seen  occuring  exactly  at  the  mid-point  of  the  conflicting 
frequencies  and  times.  These  cross-terms  build  until  they  reach  a  maximum  concentration 
at  the  transition  or  frequency  change  point.  Since  these  are  nothing  more  than  cross-terms, 
there  must  be  a  way  to  deemphasis  these  through  the  judicial  use  of  the  smoothing  window. 
Figures  13  through  16  show  the  result  of  varying  the  smoothing  window  from  a  point  of 
essentially  no  effect  to  entire  elimination  of  the  ghost  peak. 
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UNEAR  CHIRP  SINE  WAVE 
SIGNAL  LENGTH  »  0.256  SEC  /  AMP  a  1.0 


UNEAR  CHIRP  SINE  WAVE 
UNSMOOTHED  DISTRIBUTION 


Figure  12.  Unsmoothed  Wigner-Ville  Distribution 
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Figure  13.  Pseudo  Wigner-Ville  Distribution  (5  X  35  Window) 
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Figure  14.  Pseudo  Wigner-Ville  Distribution  (13  X  13  Window) 
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Amp. 


Figure  15.  Pseudo  Wigner-Ville  Distribution  (18  X  10  Window) 


Figure  16.  Pseudo  Wigner-Ville  Distribution  (35  X  5  Window) 
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This  chapter  has  shown  the  benefits  that  may  be  gained  through  the  use  of  the 
smoothing  window.  The  frequency  resolution  and  amplitude  are  the  two  major  areas  affected 
by  the  window  changes.  The  last  area  presented,  showed  the  reason  for  the  large  amplitude 
increase  at  the  frequency  transition  point  and  how  these  "ghost  peaks"  may  be  smoothed  and 
to  what  extent  deemphasized  by  adjusting  the  smoothing  window  size.  Within  the  next 
chapter,  a  thorough  energy  analysis  of  the  discretized  Pseudo  Wigner-Ville  Distribution  will 
be  conducted. 
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IV.  WIGNER-VILLE  DISTRIBUTION  ENERGY  ANALYSIS 


A.  ENERGY  DEHNITION 

As  mentioned  earlier,  the  Wigner-Ville  Distribution  is  essentially  the  Fourier  transform 
of  the  auto-correlation  of  an  input  signal.  This  may  be  thought  of  as  the  spectral  density 
function  in  a  three  dimensional  sense.  The  total  energy  of  the  Wigner-Ville  Distribution  was 
listed  earlier  as: 


where: 


1 
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j  jwDFioi ,  t)  dfjidt  - 


T|f^  -  mean  square  value  of  the  signal 


There  are  a  few  points  that  must  be  made  concerning  this  equation.  The  first  is  that  the 
integration  is  from  minus  infinity  to  plus  infinty  and  another  is  that  it  is  the  integration  of 
the  resulting  Distribution  of  a  real  signal,  not  of  an  analytic  signal  as  used  in  most  discretized 
forms  of  the  Wigner-Ville  Distribution.  Lastly  is  that  this  equation  has  been  proven 
mathematically,  which  unfortunately  does  not  take  into  account  the  fact  that  there  are 
resulting  cross-terms  and  that  the  Distribution  is  subsequently  smoothed  because  of  this. 
Since  the  last  point  is  indicative  of  what  is  seen  in  actual  real  world  problems,  the  energy  must 
be  defined  in  a  discretized  sense  that  accounts  for  the  above  mentioned  points  and  also  these 
unavoidable  cross-terms. 
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B.  EFFECT  OF  WINDOW  SIZE  ON  ENERGY 


Since  smoothing  of  the  Distribution  is  a  necessary  requirement,  the  size  of  the  applied 
window  is  something  that  must  be  decided  upon.  Another  important  point  is  knowing  what 
effect  this  will  have  on  the  energy  of  the  Distribution.  In  the  previous  chapter,  the  size 
criterion  that  was  established  by  Cartwright  [Ref.  19]  and  later  used  by  Wahl  and  Bolton  in 
Reference  12  was  discussed.  Within  this  discussion,  the  visual  benefits  of  varying  the  window 
size  in  the  time  and  frequency  domain  were  presented.  In  order  to  define  the  energy,  the 
effect  on  the  energy  of  the  Pseudo  Wigner-Ville  Distribution  as  caused  by  varying  these 
windows  will  be  presented. 

Figure  17  is  a  400  Hz  sine  wave  of  length  0.256  seconds  and  amplitude  2.0  that  will  be 
sampled  and  have  the  Wigner-Ville  Distribution  calculated.  Figure  18  is  the  resulting 
smoothed  Wigner-Ville  Distribution  that  has  been  processed  using  a  13  by  13  sized  smoothing 
window.  Figure  18  also  shows  that  the  volume  is  equal  to  0.0709118  V^.  This  is 
representative  of  the  energy  contained  within  the  Wigner-Ville  Distribution  of  the  signal  in 
Figure  17  and  is  nothing  more  than  the  volume  under  the  resulting  Distribution.  By 
definition,  this  volume  should  be  equal  to  the  mean  square  value  of  the  400  Hz  sine  wave, 
however,  due  to  the  points  mentioned  above,  this  is  not  going  to  hold  true  and  will  be 
discussed  in  depth  later.  Figures  19  and  20  are  also  smoothed  Distributions  of  this  same  400 
Hz  signal  only  using  a  35  by  5  window  in  the  case  of  Figure  19,  and  a  5  by  35  window  in 
Figure  20.  Since  N  times  M  (the  size  criterion)  are  very  close  in  these  three  examples,  the 
resulting  energy  values  of  the  three  cases  are  essentially  the  same.  The  point  to  be  emphasized 
from  this  is  that  as  long  as  the  area  of  the  window  does  not  vary  significantly,  the  energy  will 
be  represented  in  a  constant  manner  and  independent  of  the  time  and  frequency  window  size 
selected.  The  last  example  that  will  be  presented  in  this  section  is  that  this  result  is 
independent  of  the  frequency  of  the  signal  as  shown  by  the  combined  100  and  400  Hz  sine 
wave  example  of  Figure  21,  in  which  the  energy  is  equal  to  0.1418314  V^.  This  is  essentially 
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twice  the  value  of  the  single  sine  waves  of  Figures  18  through  20.  Of  special  importance  for 
transient  machinery  applications  is  Figure  22,  which  is  a  transient  or  ramped  sine  wave  that 
shows  the  same  resulting  energy  as  the  previous  stationary  signals.  This  shows  that  the 
Wigner-Ville  Distribution  energy  is  independent  of  the  nature  of  the  signal  allowing  for  the 
possibility  of  monitoring  a  piece  of  machinery  during  both  it’s  steady  state  and  transient 
phases  of  operation. 
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400  HZ  SINE  WAVE 

SIGNAL  LENGTH  =  0.256  SEC  /  AMP  =  2.0 


0.015 


400  SINE  WAVE 


Figure  19.  Pseudo  Wigner-Ville  Distribution  Energy  (35  X  5  Window) 
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400  H«  SINE  WAVE 
SMOOTHING  WINDOW  !  5  X  38 
AMPLITUDE  -  2.0 


Figure  20.  Pseudo  Wigner-Ville  Distribution  Energy  (S  X  35  Window) 


40 


0.015 


100  t  400  HZ  SINE  WAVE 
SMOOTHINQ  WINDOW  !  13  X  ia 
AMPLITUDES  -  2.0 


Figure  21.  Combined  Sine  Wave  Distribution  Energy  (13  X  13  Window) 
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Figure  22.  Ramped  Sine  Wave  Distribution  Energy  (13  X  13  Window) 
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The  previous  figures  show  that  if  the  Wigner-Ville  Distribution  is  going  to  be  used  to 
accurately  portray  or  maintain  signal  information,  then  the  window  size  must  at  least  remain 
constant  and  better  yet  to  keep  the  identical  window  size  for  each  measurement  of 
comparison.  This  requirement  to  maintain  a  constant  sized  smoothing  window  is  not  the  only 
constraint  that  must  be  meet  for  energy  analysis  purposes.  The  next  section  will  discuss 
another  of  these  elements. 

C.  TIME  DEPENDENCE  OF  WIGNER-VILLE  ENERGY 

Not  only  is  the  window  size  going  to  influence  the  energy  result,  but  the  time  length 
of  the  signal  will  also  play  a  large  part.  This  varies  considerably  from  the  conventional 
knowledge  of  the  mean  square  value  of  a  signal  being  independent  of  time.  The  standard 
mean  square  value  will  be  the  same  for  a  signal  of  time  length  2.0  seconds  as  it  will  for  a 
signal  of  4.0  seconds.  For  instance,  the  mean  square  value  of  a  sine  wave  x(t)  is; 

/v2 

(auto  spectrum)  df  -  ~ 

0 

where:  x(t)  =  X  sin  (wt) 

To  compute  the  mean  square  value  of  this  sine  wave,  it  is  independent  of  time  and  will  equal 
0.5  for  X  =  1.0.  This  is  not  the  case  for  the  mean  square  value  or  total  energy  of  the  signal 
as  found  using  the  discretized  form  of  the  Wigner-Ville  Distribution.  The  following  figures 
will  show  this  second  important  point  that  must  be  realized  when  comparing  Wigner-Ville 
energies  of  various  signals.  Figure  23  is  the  400  Hz  input  signal  of  length  0.128  seconds.  This 
time  length  is  half  of  that  of  the  400  Hz  signal  used  in  the  previous  section.  The  resulting 
Wigner-Ville  Distribution  is  shown  in  Figure  24  with  an  energy  value  of  0.0354558  V*.  As 
can  easily  be  seen,  this  resulting  energy  value  is  not  the  same  and  in  fact,  the  energy  of  Figure 
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24,  is  one  half  that  of  the  earlier  400  Hz  signals  of  length  0.256  seconds.  Had  the  signal  been 
of  a  length  0.512  seconds  or  twice  as  long  as  the  earlier  examples,  the  resulting  Wigner  Viiie 
energy  would  have  doubled.  The  window  size  used  for  this  Pseudo  Wigner-Ville  Distribution 
was  a  13  by  13  sized  window. 
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400  HZ  SINE  WAVE 

SIGNAL  LENGTH  =  0.128  SEC  /  AMP  =  2.0 


€0 

CM 


0*2  O'V  0*0  O'l.-  0*3- 


(siioA)  ganiridi^v 

Figure  23.  400  Hz  Sine  Wave  (Tl  =  0.128  sec) 
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Figure  24.  Pseudo  Wigner-Ville  Distribution  Energy  *0.128  sec) 


This  has  yielded  a  second  concern  that  must  be  kept  in  mind  when  using  the  energy  of 
the  discretized  Pseudo  Wigner-Ville  Distribution.  Since  the  energy  is  proportional  to  the 
signal  length,  an  easy  conversion  can  be  made  by  using  the  lengths  of  the  individual  raw  input 
signals.  This  time  dependence  holds  true  for  not  only  steady  state  signals  as  shown,  but  also 
for  signals  of  a  transient  nature  such  as  previously  shown  with  Figure  22.  Throughout  these 
last  two  sections,  the  comparison  of  energy  values  of  signals  has  been  mentioned  a  great  deal. 
The  thought  process  behind  this  is  that  if  there  is  an  amplitude  change  in  the  input  signal,  this 
will  be  reflected  in  the  resulting  Pseudo  Wigner-Ville  Distribution.  By  ensuring  that  these 
previous  two  conditions  are  met,  an  investigation  into  the  resulting  energy  changes  from  an 
amplitude  increase  or  decrease  of  an  input  signal  can  be  made. 

D.  ENERGY  PROPORTIONALITY  OF  THE  DISTRIBUTION 

The  mean  square  value  of  a  simple  sine  wave  is  the  amplitude  squared  divided  by  two. 
This  was  shown  in  an  earlier  equation  and  now  will  be  used  to  show  the  proportionality  of  the 
Wigner-Ville  Distribution.  For  instance,  the  following  applies  for  signals  of  differing 
amplitudes  as  shown: 

amplitude  -  1.0  -  0.5 

ampli  tude  -2.0  ijr^  -  2 . 0 

This  amplitude  increase  from  1.0  to  2.0  represents  an  increase  of  four  in  the  mean  square 
value  of  the  signal.  Figures  25  and  26  are  the  raw  input  sine  waves  of  amplitude  1.0  and  2.0 
respectively.  The  resulting  Pseudo  Wigner-Ville  Distributions  are  shown  in  Figures  27  and 
28,  which  acurately  portray  the  four-fold  energy  increase.  These  Distributions  have  been 
smoothed  using  a  window  of  13  by  13. 
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400  HZ  SINE  WAVE 

SIGNAL  LENGTH  b  0.260  SEC  /  AMP  »  1.0 


TIME  (sec) 


400  HZ  SINE  WAVE 


Figure  26.  400  Hz  Sine  Wave  (Amplitude  ■  2.0) 
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400  Hz  SINE  WAVE 
SMOOTHINQ  WINDOW  ;  13  X  13 
AMPLITUDE  =  2,0 


in 


Figure  28.  Pseudo  Wigner-Ville  Distribution  (Amplitude  *  2.0) 


SI 


These  previous  four  figures  show  the  proportionality  of  the  energy  within  the  Wigner- 
Ville  Distribution.  These  two  cases  are  just  an  example  to  show  the  proportionality  of  the 
energy.  Runs  were  made  for  many  more  cases  (ie.  X  *  4.0,  X  =  8.0,  etc)  and  the 
proportionality  relation  held  true  as  shown  with  the  previous  Figures.  This  was  another 
required  element  in  order  to  use  the  energy  as  the  basis  for  a  method  of  machinery 
monitoring.  Had  the  cross-terms  increased  in  such  a  fashion  as  to  not  allow  the  Distribution 
to  show  the  proportional  energy  change,  the  energy  content  of  the  Pseudo  Wigner-Ville 
Distribution  would  have  been  rendered  useless  for  the  purposes  of  comparison  and  evaluation. 
The  next  chapter  will  use  these  requirements  and  basics  that  have  been  developed  to  show  an 
actual  application  of  how  the  Wigner-Ville  energy  may  be  used  to  conduct  machinery 
monitoring. 
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V.  GEAR  MODEL  APPLICATION 


A.  DESCRIPTION  OF  GEAR  MODEL 

The  simple  complexity  gear  model  used  in  this  work  was  built  based  upon  a  model  used 
by  D.K.  Carlson  [Ref.  20]  for  Artificial  Neural  Network  applications  for  machinery 
monitoring.  It  may  be  of  a  simple  nature,  however,  it  provides  all  of  the  necessary 
components  that  may  be  found  in  more  complex  systems.  A  schematic  of  the  gear  model  can 
be  seen  in  Figure  29.  The  two  major  changes  of  this  model  over  the  prototype  used  by 
Carlson  is  that  the  components  are  much  larger  enabling  a  better  signal  to  noise  ratio  and 
secondly,  the  motor  used  to  drive  the  gear  model  is  mounted  on  a  separate  base  to  reduce 
motor  enduced  vibrations. 

The  model  consists  of  a  single  reduction  gear  train,  with  the  pinion  gear  being  a  90 
tooth  spur  gear  and  the  driven  gear  a  120  tooth  spur  gear.  The  gears  were  manufactured  by 
Martin  Corporation  with  a  1/2  inch  bore  and  a  14.S  degree  pressure  angle.  Aluminum  blocks 
housed  the  NICE  1/2  inch  bore  radial  ball  bearings  that  were  used  to  support  the  shaft  and 
gears.  These  aluminum  blocks  were  bolted  to  the  support  base  which  was  a  1 .0  inch  thick 
plexiglass  slab.  Driving  the  single  shaft  is  a  General  Electric  1/6^^  HP  motor,  that  is 
controlled  by  the  use  of  a  Bodine  Electric  Company  combination  rectifier  and  variable 
potentiometer  speed  controller.  The  shaft  speed  was  monitored  through  use  of  a  Power 
Instruments  Model  1 720  RPM  Optical  Proximeter.  [Ref  20] 

B.  SIGNAL  SAMPLING  EQUIPMENT 

This  section  will  describe  the  components  shown  in  Figure  30,  that  were  used  to  sample 
the  vibrational  signal  of  the  gear  model.  The  signal  was  monitored  using  an  ENDEVCO 
model  303A03  accelerometer  that  was  amplified  and  powered  with  a  PCB  model  483B07 
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power  unit.  The  signal  was  then  processed  through  a  Krohn-Hite  model  3342  analog  filter 
configured  in  a  band  pass  mode,  prior  to  being  sent  to  the  HP  3S6SS  system  for  sampling  and 
storage.  Lastly,  the  signal  was  transfered  to  the  VAX  3520  workstation  for  the  calculation  and 
plotting  of  the  Pseudo  Wigner-Ville  Distribution. 

1.  ENDEVCO  Model  303A03  Accelerometer 

The  ENDEVCO  Isotron  PE  accelerometer  is  a  minature  accelerometer  that  has  a 
medium  range  high  frequency  capability.  The  operation  of  the  accelerometer  is  based  upon 
a  piezoelectric  quartz  transducer  sensing  element.  The  sensitivity  of  the  accelerometer  is  10 
mV/g  with  a  resonant  frequency  of  70  kHz.  The  resolution  available  is  0.02  g,  with  a 
maximum  range  of  ±  500  g. 

The  accelerometer  was  mounted  on  top  of  the  aluminum  bearing  housing  that 
supports  the  shaft  and  specifically  at  the  position  shown  in  Figure  30.  The  accelerometer  wass 
attached  by  means  of  a  mounting  wax  which  yields  a  maximum  operating  range  out  to  i5k  - 
20k  Hz.  This  accelerometer  signal  was  then  sent  through  the  power  unit  described  below. 

2.  PCB  Model  483307  Power  Unit 

This  power  unit  is  used  to  power  the  low  impedance  quartz  transducers  and 
amplify  the  signal  if  desired.  The  power  unit  provides  an  adjustable  2  to  20  mA  current  for 
purposes  of  transducer  excitation.  The  gain  adjustment  is  available  from  0  to  100  and  is  set 
through  the  use  of  a  ten  turn  vernier  gain  pot  and  a  three  position  gain  multiplier  switch.  For 
the  purposes  of  these  tests,  the  gain  was  set  to  20  before  sending  the  signal  to  the  analog  filter. 

3.  Krohn-Hite  Model  3342  Anolog  Filter 

The  model  3342  variable  filter  is  a  digitally  tunned  filter  that  will  function  as  a 
High-Pass  or  Low-Pass  filter.  When  the  two  channels  of  the  filter  are  connected  together,  the 
filter  will  function  as  a  band  pass  filter,  which  is  the  configuration  for  the  gear  model  work. 
The  range  of  the  filter  is  from  0.001  Hz  to  99.9  kHz  as  adjusted  by  three  rotary  decade 
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switches  and  a  rotary  six  positiohn  multiplier  switch.  In  addition,  the  filter  unit  has  a  gain 
setting  of  unity  (0  dB)  and  10  (20  dB).  The  signal  was  not  further  amplified  using  this 
capability  of  the  filter  prior  to  sampkling  by  the  HP  3S6SS  system. 


4.  HP  3565S  Measurement  Hardware  System 

The  HP  3565S  measurement  hardware  with  the  HP-Vista  software  uses  a  HP-9000 
series  computer  system  for  controlling  purposes.  The  HP  3565S  system  is  a  modular 
multichannel  system  that  can  process  data  in  both  the  time  and  frequency  domain.  The 
system  is  capable  of  handling  64  source  and  input  modules,  but  is  presently  configured  for 
eight.  The  HP-Vista  software  allowed  for  the  sampling  and  storage  of  the  preprocessed 
accelerometer  signal  that  was  later  processed  through  the  Pseudo  Wigner-Ville 
distribution. 
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C.  ESTABLISHING  MONITORING  FREQUENCIES 


This  section  will  describe  how  the  monitoring  frequencies  have  been  determined  and 
how  they  will  be  used  to  monitor  the  gear  model  and  in  future  applications  a  piece  of 
transient  machinery.  As  shown  earlier  in  Chapter  III,  the  energy  contained  in  a  signal  is 
accurately  portrayed  within  the  resulting  time-frequency  domain  of  the  Wigner-Ville 
Distribution.  Through  the  Wigner-Ville  Distribution  energy  that  is  contained  in  the 
established  monitoring  frequencies,  faults  may  be  detected.  The  shaft  frequency  is  just  the 
speed  of  the  shaft,  while  the  gear  mesh  frequency  is  the  number  of  teeth  on  the  rotating  gear 
multiplied  by  the  shaft  speed  or  frequency.  Both  of  these  frequencies  have  associated 
harmonics  that  can  be  of  equal  importance  as  the  frequencies  themselves  in  the  n.onitoring 
of  machinery.  It  is  for  these  reasons  that  the  frequencies  listed  below  have  been  chosen  based 
upon  a  shaft  frequency  of  15  Hz. 

1.  Monitoring  Frequencies 

For  the  purposes  of  this  work,  just  the  shaft  and  gear  mesh  frequencies  of  the 
model  will  be  included.  Just  as  in  standard  established  machinery  monitoring  methods,  the 
shaft  frequency  and  associated  harmonics  can  provide  a  great  deal  of  information  for  fault 
detection  and  diagnosis.  Along  with  the  shaft  frequencies,  the  gear  meshing  frequencies 
provide  valuable  information  of  faults  and  will  also  be  monitored.  The  purpose  of  the  two 
examples  will  be  to  show  that  the  Wigner-Ville  Distribution  energy  (ie.  volume)  changes  as 
a  result  of  damage  and  will  prove  to  be  a  very  valuable  tool  in  future  machinery  monitoring 
methods.  The  established  monitoring  frequencies  have  been  filtered  as  follows: 

•  Shaft  Input  1:  (S  -  100  Hz)  This  is  the  spectrum  of  signals  in  the  lower  frequency  range 
selected  to  determine  the  total  energy  contained  within  the  shaft  frequency  and  it’s 
harmonics. 

•  Shaft  Input  2:  (10  -  20  Hz)  This  frequency  range  captures  just  the  shaft  frequency. 

•  Shaft  Input  3:  (25  -  35  Hz)  This  frequency  range  captures  the  1*‘  shaft  harmonic. 
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Shaft  Input  4:  (40  -  50  Hz)  This  frequency  range  captures  the  2"'^  shaft  harmonic. 


•  Shaft  Input  5:  (55  -  100  Hz)  This  frequency  range  captures  the  upper  harmonics  of  the 
shaft  frequency. 

•  Gear  Mesh  Input  1:  (1250  -  1450  Hz)  This  frequency  range  captures  the  gear  meshing 
frequency. 

•  Gear  Mesh  Input  2:  (2600  -  2800  Hz)  This  frequency  range  captures  the  1**  gear  mesh 
harmonic. 

•  Gear  Mesh  Input  3:  (3950  -  4150  Hz)  This  frequency  range  captures  the  2"*^  gear  mesh 
harmonic. 


2.  Combined  Monitoring  Technique 

Using  the  Wigner-Ville  energy  as  an  input  to  an  Artificial  Neural  Network  is  a 
follow-on  goal  of  this  research.  It  is  for  this  reason,  that  the  above  established  frequencies 
have  been  identified  as  inputs.  Within  this  scheme,  the  Wigner-Vilie  is  used  to  characterize 
the  signal  and  give  a  single  energy  value  at  the  particular  frequency  of  interest  as  an  input  to 
an  Artificial  Neural  Network.  By  filtering  these  signals  prior  to  processing,  just  the 
frequency  of  interest  will  be  captured  along  with  it’s  respective  Wigner-Ville  energy  content. 
The  Wigner-Ville  energy  data  can  now  be  used  to  produce  a  comprehensive  training  set  for 
the  Artificial  Neural  Network.  By  pre-processing  the  inputs  to  a  Neural  Network  with  the 
filtering  and  the  Wigner-Ville  Distribution,  the  number  of  inputs  is  greatly  reduced.  Since 
these  associated  energy  levels  are  fairly  stable  for  a  given  machinery  condition,  the  Neural 
Network  will  be  capable  of  detecting  patterns  in  energy  level  shifts  that  are  associated  with 
particular  machinery  faults.  To  support  this  last  statement,  the  following  section  will  show 
two  examples  of  machinery  faults  imposed  on  the  gear  model  and  the  response  of  the  Wigner- 
Ville  energy  to  these  faults. 
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D.  MACHINERY  FAULT  EXAMPLES 
1.  Mechanical  Looseness 

This  fault  was  imposed  on  the  gear  model  by  loosening  the  aluminum  support 
blocks  that  house  the  bearings  and  support  the  shaft.  Specifically,  the  blocks  that  support  the 
driver  shaft  were  loosened.  The  degree  of  looseness  imposed  was  difficult  to  judge  and  made 
this  particular  fault  a  difficult  one  to  diagnose.  The  blocks  were  not  completely  undone,  but 
loosened  just  a  small  amount  from  their  original  firmly  bolted  posture.  Shown  in  the  table 
below  are  the  results  of  averaging  two  good  runs  and  two  post-damage  runs  and  comparing 
the  respective  energy  values  at  the  different  frequencies.  These  energy  values  are  nothing 
more  than  the  volume  under  the  Pseudo  Wigner-Ville  Distribution  as  discussed  earlier.  Only 
the  shaft  frequencies  are  shown  for  purposes  of  clarity  and  since  it  is  in  these  monitored 
frequencies  that  a  fault  of  mechanical  looseness  should  present  itself. 
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T*ble  I.  MECHANICAL  LOOSENESS 


MECHANICAL  LOOSENESS 

ENERGY  VALUES  (V"2) 


Freouencv 
Shaft  Spectrum 
Shaft  Frequency 
1st  Harmonic 
2nd  Harmonic 
Upper  Spectrum 


Good 

Damaoed 

41.0 

28.0 

1.07 

0.83 

2.31 

1.22 

1.23 

2.40 

31.4 

26.31 

As  can  be  seen  in  the  table  above,  the  energy  content  approximately  doubled  at  the  second 
shaft  harmonic.  This  is  a  very  encouraging  re:>ult  for  the  purposes  of  using  the  Pseudo 
Wigner-Ville  Distribution  energy  as  a  fault  detection  tool. 

2.  Gear  Tooth  Damage 

This  particular  fault  was  imposed  on  the  gear  model  by  shaving  a  gear  tooth  on 
the  driver  gear  The  degree  of  damage  was  much  easier  to  predict  in  this  case  and  made  this 
particular  fault  an  easier  one  to  control.  The  gear  tooth  that  was  damaged  had  approximately, 
l/8‘^  of  the  leading  edge  removed.  Shown  in  the  table  below  are  again  the  results  of 
averaging  two  good  runs  and  two  post-damage  runs  and  comparing  the  respective  energy 
values  at  the  different  frequencies.  For  this  case,  just  the  frequencies  associated  with  the  gear 
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mesh  are  presented,  since  it  is  in  these  monitored  frequencies  that  a  fault  involving  the 
damage  of  a  gear  train  would  present  itself. 

Table  II.  GEAR  TOOTH  DAMAGE 


GEAR  TOOTH  DAMAGE 


ENERGY  VALUES 

< 

> 

MESH  FREQUENCY 

GOOD 

DAMAGED 

Fundamental 

17.14 

154.4 

1st  Harmonic 

0.53 

3.44 

2nd  Harmonic 

0.16 

0.48 

As  shown  in  the  above  table,  the  energy  content  greatly  increased  at  the  fundamental  gear 
mesh  frequency.  Additionaly,  there  was  also  increases  at  the  harmonics  and  shown  in  Figures 
3 1  and  32  is  a  graphical  presentation  of  the  energy  increase  that  occured  due  to  the  damaged 
tooth. 
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Figure  31.  Gear  Mesh  Frequency  (Pre-Damage) 
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This  chapter  has  shown  two  examples  of  machinery  damage  and  the  resulting  changes 
of  the  Pseudo  Wigner-Ville  Distribution  energy.  By  applying  the  window  and  energy 
requirements  discussed  in  Chapters  III  and  IV,  the  Wigner-Ville  energy  has  yielded  a  definite 
change  to  allow  a  means  for  machinery  monitoring.  The  Wigner-Ville  energy  will  make  a 
stable  input  to  a  Neural  Network  and  provide  a  quick  efficient  method  by  pre-processing  the 
machinery  signal  with  the  Wigner-Ville  Distribution  as  shown.  This  pre-processing  will  allow 
the  number  of  inputs  to  the  Neural  Network  to  remain  at  a  reasonable  level  and  provide  for 
a  manageable,  quicker  hardware  system. 
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VI.  CONCLUSIONS 


The  Pseudo  Wigner-Ville  Distribution  has  been  shown  to  provide  an  excellent 
characterization  of  an  input  signal.  Based  upon  the  research  conducted,  it  will  also  provide 
a  means  for  machinery  condition  monitoring  and  diagnostics.  The  research  has  also  shown 
that  in  order  to  use  the  Wigner-Ville  Distribution  in  a  machinery  monitoring  scheme,  the 
following  conditions  or  conclusions  must  be  met  and  understood. 

•  The  smoothing  window  size  applied  must  remain  constant  for  purposes  of  comparing 
various  signals. 

•  The  individual  time  and  frequency  lengths  of  the  smoothing  window  may  be  adjusted 
to  provide  a  better  frequency  resolution  or  constant  amplitude. 

•  The  energy  contained  within  the  Pseudo  Wigner-Ville  Distribution  of  a  signal  is 
proportionaly  time  dependent. 

•  The  energy  of  the  Distribution  is  proportional  to  the  mean  square  value  of  the  input 
signal  and  is  not  adversely  affected  by  the  cross-terms. 

•  L?  tly,  the  damaged  gear  model  conditions  showed  a  definite  change  in  the  energy  of 
the  Pseudo  Wigner-Ville  Distribution  at  the  respective  monitoring  frequencies. 
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VII.  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 


The  Pseudo  Wigner-Ville  Distribution  computer  program  that  was  established  here  at 
the  Naval  Postgraduate  School  is  very  stable.  The  pronounced  peaks  at  the  transition  points 
are  now  definable  and  additionally,  with  the  energy  of  the  Pseudo  Wigner-Ville  Distribution 
now  being  a  definable  and  understood  quantity,  the  next  additional  research  in  this  area  may 
include: 

•  Establish  a  baseline  or  data  bank  of  "good"  Wigner-Ville  Distribution  data  for  the 
present  gear  model. 

•  Conduct  additional  runs  of  various  levels  of  damage  of  gear  model  components  for 
comparison  with  good  data. 

•  Establish  a  Torpedo  Ejection  Pump  baseline  of  good  data  as  with  the  gear  model. 

•  Continue  with  the  research  of  linking  the  Pseudo  Wigner-Ville  Distribution  as  an  input 
to  a  Neural  Network  Application. 
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APPENDIX  COMPUTER  PROGRAMS 

The  following  programs  have  all  been  developed  here  at  the  Naval  Postgraduate  School 
and  used  in  the  conduct  of  the  research  presented.  The  language  used  is  FORTRAN  77  with 
the  plotting  programs  using  calls  to  CA-DISSPLA  software. 


A.  WVD3  COMPUTER  PROGRAM 

PROGRAM  UVD3 

THE  FOLLOWING  INDIVIDUALS  HAVE  BEEN  INSTRUMENTAL  IN  THE 
DEVELOPMENT  OF  THIS  PROGRAM  AND  PREVIOUS  VERSIONS: 

Graham  Rossano 
Jose  G.  Calahorrano 
Scott  G.  Spooner 
Prof.  Y.S.  Shin 


This  program  calculates  the  Wigner  Distribution 
Function  of  a  signal.  It  includes  the  option  of  applying  a 
smoothing  function  and  reducing  the  distribution  to  different 
sizes  while  writing  the  results  to  files  for  subsequent 
plotting. 


VERSION  3 

VARIABLE  DECLARATION 

************* 

-  setting  nuiR)er  of  data  points  ■■ 


integer  43,dp2,mvopt,redopt,dimt,dimf ,nn,nn,rf ,rdp,rdp2, 

:  nf,mt,nf2,mt2 

PARAMETER  (DP  ■  S12) 

real  pi ,tin(dp},ain(dp),rawa(dp),rwdf(dp*2,dp),r8wdf(dp*2,dp}, 
r  awt  ( dp) ,  dt  sun,  de  1 1 ,  dt ,  asuRi,  meanv ,  mt  i  me ,  de  1 1 ,  de  1 2 , 
const,  t,  sun,  SkMb,val,sval,wdf(dp*2,Cp),coef,df,  AREA 
coaplex  s(dp*2),dun,c(dp*5),duii3,dun2 
character*25  inname 

common/ in/  iHvopt,redopt,dimt,dimf ,inn,nn,rf ,rdp,rdp2. 
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nf,Mt,nf2,iit2 

coMMn/rl/  pi,tin,ain,rMM,rwdf ,r>udf , 

r  awt ,  dt  tuR ,  de  1 1 .  dt ,  asun,  iManv ,  lit  1 M ,  de  11 ,  de  1 2 , 
const ,  t ,  Sin,  ti«b,  va  I ,  sva  I ,  wdf ,  coef  ,df ,  area 
coamon/coap/  s,dua,c,dkn3,duii2 
dfa  «  dp*2 


INITIALIZING  VARIABLES 


—  Print  description  of  prograai  — 
print* 

print*,'  Prograai  to  calculate  the' 
print*,'  Pseudo  UIgner-Ville  Distribution' 
print* 

---  Set  initial  values 

print*,'  Enter  nasK  of  signal  input  file* 

read(5,90A)  InnaaK 

print* 

print*,'  Do  you  wish  to  reaiove  the  aiean  value?' 
print*,'  Enter  1  for  yes  -  0  for  no* 
read(5,910)  aivopt 
print* 


print*,'  Input  the  reduction  size  desired' 
print*,*  Input  1  for  64  by  32  * 


print*,'  input  2  for 
print*,'  input  3  for 
print*,'  Input  4  for 
read(S,910)  redopt 
print* 


128  by  64 
128  by  128 
256  by  128 


print*,'  Input  the  smoothing  window  size  desired* 
print*,'  Input  the  size  for  the  frequency  axis' 
read(5,910)  nf 

print*,'  Input  the  size  for  the  time  axis' 
read(5,910)  mt 


calculation  of  some  coamon  used  constants 


pi  *  4.0  *  atan(I.O) 


*  THIS  IS  THE  CALCULATION  PART  OF  THE  PROGRAM 


open(4,f I le>Inname,status>'old' ) 
rewind  4 


*  PRELIMINARY  CALCULATIONS: 


CALL  0ATAIN3 


wrlte(6,*)  ,  '  done  with  datalrwS* 


CALL  0TCALC3 
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write(6,*)  ,  '  done  with  dtcalc3' 

CALL  MEAN3 

Hrite(6,*)  ,  '  done  with  iiean3' 
write(6,*)  ,  '  Mean  value  reaioved  Meanv 

*  SIGNAL  MODIFICATIONS: 

*  Window  application  (Modified  haaning  window) 

CALL  HAMMG3 

write(6,*)  ,  '  done  with  haamg3' 

*  Converaion  of  real  signal  to  an  analytical  one 

CALL  ANAL3 

write(6,*)  ,  ■  done  with  anal3' 

*  CALCULATION  OF  THE  UIGNER  VILLE  DISTRIBUTION: 

CALL  WIGH3 

write(6,*)  ,  '  done  with  wigh3' 

*  REDUCTION  t  SMOOTHING  OF  THE  WIGNER  VILLE  DISTRIBUTION: 

CALL  REDUCE3 

write(6,*)  ,  '  done  with  reduce3' 

CALL  SMOOT H3 

t/rite(6,*)  ,  '  done  with  sinooth3' 


WRITING  OF  THE  DIFFERENT  ARRAYS  GENERATED 
TO  FILES  FOR  PLOTTING  PURPOSES. 


OPENING  OF  THE  OUTPUT  FILES 

OPEN(UH I T*7, FILE*' RAWT I ME .OUT ■ , STATUS* 'NEW' ) 
OPEN(UN I T *B . F I LE* ' WNDWT I M . OUT ' , STATUS* ' NEW ' ) 
OPEN(UN I T*9, F I LE* ' RWDF .OUT ' , STATUS* ' NEW' ) 
OPEN( UN I T* 1 0 . F I LE* ' RSUD F . OUT ' , STATUS* ' NEW ' ) 
OPEN(UNIT*11,FILE*'WDF.OUT',STATUS='NEW') 


•  WRITING  OF  RAW  AND  WINDOWED  TIME  SIGNAL  TO  FILES 

DO  300  I  *  I, OP 
WRITE(7,930)  TIN(I),RAWA(I) 

WRITE(8,930}  TIN(I),AIN(I) 

300  CONTINUE 

•  WRITING  OF  REDUCED/UNSMOOTHED  WVD  TO  FILE 

DO  400  I  *  1, dp/mm 
DO  400  J  «  1,DP2/nn 
WRITE(9,906)  RWDF(J,I) 

400  CONTINUE 

•  WRITING  OF  REDUCED  t  SMOOTHED  WVD  TO  FILE 
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* 

DO  500  I  X  1,DP/Mii 
DO  500  J  «  1 ,DP2/nn 
URITE(10,906}  RSUDF(J,1) 
500  CONTINUE 

•  WRITING  OF  ENTIRE  UDF  TO  A  FILE 

* 

DO  600  I  «  1,DP 
DO  600  J  >  1 ,DP2 
URITE(11.906)  UDF(J,I) 
600  CONTINUE 

*  FortMt  ttatanents 

904  fortMt(a2S) 

906  fonnat(2X,g16.8) 

910  formt(i6) 

930  fonMt(2x,g16.8,5x, 016.8} 

END 

*  SUBROUTINES 

include  'ANALS.INC' 
include  'CORRB.INC 
include  'DATAINS. INC 
include  'DTCALCB.  INC 
include  'FFTB.INC 
include  ‘HAMMGS.INC 
include  'MEANS. INC 
include  'SMOOTHB.INC 
include  'REDUCES. INC 
include  'UIGHS.INC 
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B.  SUBROUTINES  OF  WVD3 


Listed  below  in  alphabetic  order  are  the  subroutines  used  in  the  main  program  WVD3. 


SUBROUTINE  ANAL3 

r 

>*********************«*«*****************************«**1 


This  subroutine  converts  a  real  signal  to  an  ' 

analytic  one  ’ 

1 

r*******«*************«**********«****««***«*«***«*#***i 


integer  dp.dpZ.mvopt.redopt.dimt.dimf .itn.nn.rf ,rdp,rclp2, 

:  nf  ,nit,nf2,mt2 

PARAMETER  (OP  >  512) 

real  pi ,tin(dp),ain(dp),rawa(dp),rudf (dp*2,dp),rsudf(dp*2,dp), 
rawt  (dp) ,  dt  sun,  de  1 1 ,  dt ,  asun.  meanv,  mt  i  me ,  de  1 1 ,  de  1 2 . 
cons t . t , sun, sunb, va I , s va I . wdf (dp*2 , dp) , coef , df , area 
complex  s(dp*2),dun,c(dp*5),dun3,dun2 
character*25  inname 

common/ in/  mvopt,redopt,dimt,dimf ,nin,nn,rf ,rdp,rdp2, 

:  nf ,mt,nf2,mt2 

common/rl/  pi,tin,ain,raua,rwdf ,rswdf , 

r  awt ,  dt  sun,  de  1 1 ,  dt ,  asun,  meanv ,  mt  i  me ,  de  1 1 ,  de  1 2 , 

:  can8t,t,sun,sunb,val,sval,wdf ,coef ,df,area 

connon/camp/  s,dun,c,dun3,dun2 
dp2  >  dp*2 


do  100  i>1,dp 
suribO.O 

do  200  j>1,dp 
aui^O.O 

if(i-j.eq.O)  go  to  200 
n«i-j 

val*pi*n/2. 

sval>sin(val) 

sunb^a i n( j )*sva I *sva I / va I 

200  susBSUtH’Suit) 

s(i)«cmplx(ain(i),sun) 

100  continue 


END 
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SUBROUTINE  DATAIN3 

•A*********************************************************** 

* 

This  subroutine  assuses  the  input  file  is  in  the  following  * 
format:  time  amplitude  (2x,e16.8,5x,e16.8)  with  an  end  of  * 
file  indicator  of  a  last  entry  9999., 9999.  * 

* 

**«*****#♦*«******«****«•*♦*•***«••**♦************•******«•** 


integer  dp,dp2,myopt,redopt,dimt,dimf ,mm,nn,rf .rc^.ro^, 

:  nf  ,ait,nf2,mt2 

PARAMETER  (OP  «  512) 

real  pi,tin(dp),ain(dp),rawa(dp).rwdf(dp*2,dp),rswdf(dp*2.dp), 
r  awt  ( dp) ,  dt  sun,  de  1 1 ,  dt ,  asutt,  meanv ,  ait  i  me ,  de  11 ,  de  1 2 , 
cons t , t , sum, sunb, va I , 8 va I , wdf ( dp*2 , dp) , coef ,df , ARE A 
complex  s(dp*2),dun,c((^*5),di^,dun2 
character*2S  inname 

common/in/  mvopt,redopt,dimt,dimf ,mm,nn,rf,r(^,rdp2, 
nf ,mt,nf2,mt2 

conmon/rl/  pi ,tin,ain,rawa,rwdf ,rswdf , 

r  awt ,  dt  sun ,  de  1 1 ,  dt ,  asun,  meanv ,  mt  i  aie ,  de  1 1 ,  de  1 2 , 
const ,  t ,  sun,  sunb,  va  I ,  sva  I ,  wdf ,  coef ,  df  ,area 
conmon/comp/  s,dun,c,duii3,dun2 
dp2  >  dp*2 


'simple  loop  to  read  in  time  8  amplitude’ 
do  100  j>1,dp 

re8d(4,902)  rawt(j)  ,  rawa(j) 
tin(j)  *  rawt(j) 
ain(j)  *  rawa(j) 


100  continue 

*  FORMAT  STATEMENT 

902  FORMAT(2X,E16.8,SX,E16.8) 

200  END 
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SUBROUTINE  DTCALC3 


>•**••••••****•*•*•«•••••**•**•••*••**•******••*••••• 

This  subroutine  calculates  the  delta  t  of  the  signal 


integer  clp.clp2,iirvopt,redopt,diiitt,dinif ,iin,nn,rf,rdp.rdp2, 

:  nf ,mt,nf2,mt2 

PARAMETER  (DP  =  512) 

real  pi,tin(dp),ain(dp),raua(dp)>rwdf(dp*2.dp),rswdf(dp*2,dp), 
raHt(dp),dtsijn,delt,dt,asun,iaeanv.mtime,del1.del2, 
cons t , t . sun. sumb, va I , sva I , wdf (dp*2 , dp) , coef , df , ARE A 
complex  s(dp*2),dun,c(dp*5),dun3,dun2 
character*25  inname 

connon/in/  mvopt.redopt.dimt.dimf .mn.nn.rf .rdp.rdp2. 

:  nf ,mt,nf2,mt2 

cammon/rl/  pi.tin.ain.raua.rwdf .rswdf , 

;  rawt,dtsun,delt,dt,asun,meanv,intiine,del1,del2, 

:  const.t.sum.sunb.val.sval.wdf ,coef ,df .area 

comnon/camp/  s.dijn.c,dutn3.duni2 
dp2  =  dp*2 


* 


dtsun  s  0.0 

do  100  i  =  1  ,  dp-1 

delt  =  tin(i+1)  -  tin(i) 
dtsun  *  dtsun  ♦delt 

1 00  cont i nue 

dt  *  dtsun  /  (dp-1.) 

END 
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SUBROUTINE  FFT3 


This  subroutine  calculates  the  Fast  Fourier  Transform  * 

* 


integer  <^,dp2,nivopt,redopt,difflt,dimf .nHjnn.rf .rdp,rdp2, 
nf ,mt,nf2,mt2 
PARAMETER  (DP  «  512) 

real  pi,tin(dp),ain(dp).raua(dp).rudf(dp*2,c^>,rswdf(dp*2.dp), 
r  aw  t  ( dp ) ,  dt  sun,  de  1 1 ,  dt ,  asm,  Rieanv .  mt  i  sie,  de  11 .  de  1 2 , 
const ,  t ,  sm,  smb,  va  I ,  s  va  I ,  wdf  (dp*2  ,dp) ,  coef  ,df 
complex  s(dp*2),dm,c(dp*5),dm3,dm2 
character*25  inname 

cannon/ in/  invopt,redopt,difflt,dimf ,inn,nn,rf,rdp,rdp2, 

:  nf ,mt,nf2,mt2 

common/rl/  pi ,tin,ain,raua,rwdf ,rsudf , 

rawt,  dtsm,  del  t,dt,  asm,  meanv,mtiine, dell, del2, 

:  const ,  t ,  sm,  smb,  va  I ,  sva  I ,  udf ,  coef  ,df ,  area 

cornnon/comp/  s,dm,c,dun3,dutn2 
dp2  >  dp*2.0 


constsdp2 

val3alog(const)/alog(2.)>.1 

j=1 

JO  40  i»1,dp2-1 
if  (i.ge.j)  9o  to  10 

dm3=c(  j) 
c(j)*c<i) 
c(i  )=dun3 
10  k=dp 

20  if  (k.ge.j)  go  to  30 

j=j-k 
k=k/2 
go  to  20 
30  j-j*k 

40  cont i nue 

do  70  n«1,val 
coef=2**n 
coef=coef/2 
duni2=cfflplx(1.,0.} 

duDBcmplxC cos( pi /coef ),- s i n(pi /coef ) ) 

do  60  j*1,coef 
do  50  i*j,dp2,2*coef 
i i^i+coef 
dm3=c(  i  i  )*dm2 
c{ii)=c(i)-dum3 
c{  i  >=c( !  )*dum3 
50  cont i nue 

dm2»dm2*dm 
60  continue 
70  continue 
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SUBROUTINE  NAMMG3 


* 


*  This  subroutine  applies  a  siodified  hanning  windou 

*  to  the  >.)gnal  ain(t> 


* 

integer  4>,db2,iiivopt,redopt,diint,diarf,iiin,nn,rf,r4>,re^, 

:  nf ,Mt,nf2,mt2 

PARAMETER  (DP  «  512} 

real  pi . t i n(dp) , a i n(dp) , rawa(dp) , rudf (dp*2,dp) , rswdf (dp*2,dp) , 
raut(dp),dtsiin,delt,dt,asijn,nieanv,iiitime,del1,del2, 
cons t , t , sun, sumb, va I , sva I , wdf ( dp*2 ,dp) , coef ,df , area 
complex  s(dp*2>,dun,c(dp*5),duR3,duRi2 
character*25  irmame 

* 

comnon/in/  nivopt.redopt.difflt.dimf ,iiin,nn,rf ,rdp,rdp2, 

:  nf  ,int,nf2,mt2 

connon/rl/  pi,tin,ain,rawa,rwdf ,rswdf, 

:  raut,dtsun,delt,dt,asun,ineanv,mtiine,del1,del2, 

:  const ,  t ,  sun,  sunb,  va  I ,  sva  I ,  udf ,  coef  ,df ,  area 

cornnon/comp/  s,dun,c,dum5,dun2 
dp2  >  dp*2 

* 


mtiniesdp*dt 

del1=0.1*nitime 

del2=0.9*mtiine 

const=pi/del1 

do  100  j  <  1  ,  dp 
t  «  <j-1)  *  dt 
if  (t.le.dell)  then 

ain(j)  =  8in(j)  •  < .5*(1 .-cos(const*t))) 
elseif  (<t.ge.del2).and.(t.le.intime})  then 
ain(  j)xain(  j)*(  .5*(1  .-cos<const*<intime-t>>>) 
else 
endif 

100  continue 
END 
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SUBROUTINE  NEAN3 


'  This  subroutine  calculates  and  removes  the  mean  value  ’ 
'  of  the  signal.  ' 

r  ' 

r********************************************************' 


integer  dp,dp2,mvopt.redopt,dimt,dimf ,mn.nn,rf,rdp,rdp2, 
nf ,mt,nf2,mt2 
PARAMETER  (OP  >  512) 

real  pi,tin«^>),ain(dp),rai(a(dp),rtidf<4)*2,c^»),rsMdf(dp*2,dp), 
rawt(dp)  .dtsun.del  t  ,dt  .asun.meanv.mt  ime.dell  ,del2, 
cons t , t , SUP, suMb, va I , s va I , wdf (dp*2 .  ) , coef , df , AREA 

complex  B(dp*2),din,c(d^*S),dum3,diJii2 
character*2S  inname 

* 

common/in/  mvopt,redopt,dimt,dimf ,iMi,nn,rf,rdp,r4)2, 
nf ,mt,nf2,mt2 

conmon/rl/  pi,tin,ain,raMa,rwdf,rsudf, 

r  aut ,  dt  sum,  de  1 1 ,  dt ,  asm,  meanv,  mt  i  me ,  de  1 1 ,  de  1 2, 

:  const,  t,sm,sumb,val,sval,udf,coef,df,  area 

cooiDon/comp/  s,dm,c,dm3,dm2 

dp2  X  dp*2 

* 

* 

asm  =  0.0 

do  100  i  X  I  ^  dp 

asm  X  asm  *  ain(i} 

100  continue 

meanv  x  asm  /  dp 

if  (mvopt.eq.l)  then 

do  200  i  X  1  ,  dp 
oin(i)  X  ain<i)  •  meanv 
200  cont i nue 

endif 

END 
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SUBROUTINE  REDUCES 


This  subroutine  reduces  the  wigner  viUe  ' 

to  a  workable  size. 


integer  elp,dp2,mvopt,redopt,dimt,diiBf ,ani,nn,rf,rdp,rdp2, 

:  nf  ,int,nf2,nt2 

PARAMETER  (DP  «  512) 

real  pi , t in(dp) , ain(dp) .  rawa(dp) , rwdf (dp*2,dp) , rswdf (dp*2,^> , 
rawt(dp),dtsun,delt,dt.asiin,iiieanv,MtiMe,del1,del2j 
const , t , sua, suib, va I , sva I , wdf (dp*2 ,dp) ,coef ,df , area 
coaplex  s(dp*2),dun,c(dp*5),dua3,dijii2 
character*25  imame 

conmon/in/  aivopt,redopt,di«t,di«if ,aMi,nn,rf ,r^,rdp2, 
nf  ,mt,nf2,iiit2 

conmon/rl/  pi,tin,ain,rawa,rwdf,rswdf , 

rawt , dt sum, de I t , dt , asum, meanv, mt i me. de I 1 , del 2 , 
const , t , sun, sumb, va I , sva I , wdf , coef , df , area 
conmon/cotnp/  s,dun,c,dun3,dum2 
dp2  >  dp*2 


IF  (REDOPT.EQ.1)  THEN 
RF  «  64 
RT  »  32 

ELSEIF  (RED0PT.E0.2)  THEN 
RF  «  128 
RT  «  64 

ELSEIF  (RED0PT.E0.3)  THEN 
RF  «  128 
RT  >  128 
ELSE 

RF  =  256 
RT  =  128 

END  IF 

rm  -  dp2  /  RF 
mm  z  dp  /  RT 

I  *  0 

do  100  j  =  1  ,  dp2  ,  nn 
1*1  +  1 
k  *  0 

do  100  i  *  1  ,  cl^  ,  mn 
k  *  k  +  1 

rwdfd.k)  »  wdf<j,i) 
100  continue 

END 
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SUBROUTINE  SMOOTHS 


'  This  subroutine  reduces  and  smooths  the  UDF  along  both 
'  the  time  and  frequency  axes  for  purposes  of  plotting 
'  clarity. 

r 

r*********************************************************** 

t 

real  hg(-100:100, -100:100) 


integer  dp,dp2,mvopt,redopt,dimt,dimf ,iiis,nn,rf,rdp.rdp2, 
nf  ,sit,nf2,mt2 
PARAMETER  (DP  >  512) 

real  pi ,tin(dp),ain(dp),rawa(dp),rudf ((^*2,o^),rsudf(dp*2,dp), 
r awt ( dp) , dt sun, de 1 1 , dt , asum,  meanv ,  mt  i  me, de 1 1 . de 1 2 . 
cons  t ,  t ,  sum,  sunb ,  va  I ,  s  va  I ,  udf  ( e^2 ,  dp  ) ,  coef ,  df ,  a  rea 
complex  8(dp*2),dun,c(dp*5),di^3,dijn2 
character*2S  inname 

common/  i  n/  mvopt ,  redopt ,  di  mt ,  di  mf  ,mm,  m,  rf ,  rdp,  rdp2, 

:  nf ,mt,nf2,mt2 

common/r I /  pi , t i n, ai n, rawa, rudf , rsudf , 

rawt , dt sun , de 1 1 , dt , asun,  meanv , mt i me , de 1 1 , del 2 , 

:  const, t, sun, 8umb,val,8val,wdf,coef,df, area 

comnon/cofflp/  s,dun,c,dun3,dum2 
dp2  =  dp*2 


df=1./<A.*dp*dt) 

rdpodp/nm 

rdp2sdp2/nn 

nf2«nf*2 

mt2»mt*2 

val=1 ./(<2.*pi )**2*nf*mt) 

do  20  j«-mt2,mt2 
do  10  i»-nf2,nf2 

coef*-((j*j)/<2.*mt*mt)+<i*i)/(2.*nf»nf)) 
hg( i , j )*val*exp(coef ) 

10  continue 
20  continue 

do  100  j>1,rdp 
do  100  i«1,rdp2 
100  rswdffi, j)>0.0 


ii  ■  0 

DO  500  N3l,DP2,nn 
jj  »  0 
ii  »  ii*1 

DO  450  M  X  1,DP,inm 
jj  “  jj+1 

DO  400  K>m-MT,m«MT 

DO  350  L  *  n-nf,n*nf 

IF  (L.LT.1)  THEM 

RSVBF(ii,jj)  «  RSW0F<ii,jj)  ♦  0 
ELSEIF  (K.LT.1)  THEN 
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RSUDFCii.jj)  »  RSWDF(ii.jj)  ♦  0 
ELSEIF  (L.GT.DP2)  THEN 

RSWDFdi.jj)  *  RS«)F(ii,jj)  ♦  0 
ELSEIF  (K.GT.DP)  THEN 
RSWDFdi.jj)  «  RSWDFdiJj)  ♦  0 
ELSE 

RSWDFdi.jj)  -  HSWDFdi.jj)*«)F(L.K)*HG{L-N,K-M) 
ENOIF 

350  CONTINUE 

400  CONTINUE 

RSUDF( I f . i i )sRSUDF< ii.jj) 

450  CONTINUE 

500  CONTINUE 

END 
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SUBROUTINE  UIGH3 


t*****************************************i 


This  subroutine  calculates  the  UDF  of  the  signal  * 

* 

r**************************************************** 


integer  dp,dp2,iavopt,redopt,dint,disif ,mi,nn,rf,rdp.rdp2, 

:  nf,iiit,nf2,flit2 

PARAMETER  (DP  >  512) 

real  pi,tin(dp),ain(d^),raua(dp),rwdf(dp*2,dp),rsMdf(dp*2,dp>, 
r  awt  ( dp) ,  dt  suR.  de  1 1 ,  dt ,  asum,  meanv .  lat  i  me,  de  11 ,  de  1 2 , 
const ,  t ,  sun.  SMb,  va  I ,  sva  I ,  Mdf  (dp*2,  dp) ,  coef  ,df ,  area 
complex  s(dp*2).diM,c((^*5),dun3.dijin2 
character*2S  imame 

comnon/in/  mvopt.redopt.dimt.dimf ,Rm,nn.rf,rdb.rdp2, 
nf,fflt,nf2,iRt2 

common/rl/  pi.tin.ain.raua.rwdf.rsudf, 

r  aut ,  dt  sum.  de  1 1 ,  dt ,  asm.  meanv .  mt  i  me .  de  1 1 ,  de  1 2 , 
const .  t .  sm.  srnb.  va  I .  sva  I .  wdf .  coef  .df .  area 
common/comp/  s.dm.c.dijii3.dm2 
dp2  -  dp*2 


do  100  j  3  I  ,  dp 

*******auto  correlation  of  the  signal  corr3 
coef  *  2.0  *  dt 
do  300  i  >  1  .  dp»1 

if(j.ge.i)  then 
dm  «  s(j-i>1) 
else 

dm  «  cmplx(0..0.) 
endif 

c(i)  ■  coef  *  <s(j*i-1)*conjg(dm)) 

if(i.ne.j.afKl.i.ne.dp»1)  then 
c(dp2-i'»2)  »  conjg(c(i)) 
endif 

300  continue 

call  fft3 

do  200  i>1.4>2 

*idf(i.  j)»real(c(i)) 

200  continue 

100  continue 
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C.  PLOTTING  PROGRAMS 


Listed  below  are  the  two  computer  programs  that  were  used  to  produce  the  Figures  in 


this  work. 

1.  2-D  Plotting  Routine 


PROGRAM  SIGNALPLT 

•A******************************************************* 

• 

*  THIS  PROGRAM  USES  THE  GRAPHIC  PACKAGE  DISSPLA  TO 

*  PLOT  THE  A  2-D  FUNCTION. 

* 


•  ***-- -declaring  variables---*** 

REAL  XARAY(1000),YARAY(1000) 

RE AL  X . Y . F , XST ART , XEND , YSTART , YENO , XSTEP 
INTEGER  IPROPT.NPNTS.IX.IY 

* 

*  ***---INITIALI2ING  VARIABLES  AND  SETTING  GRAPH  OPTIONS---*** 

* 

*  — set  option  for  x  starting  &  ending  value — 
WRITEC*,*). 'DECIDE  THE  STARTING  AND  ENDING  VALUES  OF  X' 
URITE(*,*),' INPUT  THE  STARTING  VALUE  OF  X* 

READ(*,*),  XSTART 

WRITEt*,*), 'INPUT  THE  EHDiHG  VALUE  OF  X' 

READ!*,*),  XEND 

*  ---set  desired  x  step  siie--- 
URITE(*,*), 'DECIDE  THE  X  STEP  SIZE' 

UR1TE(*,*),' INPUT  THE  X  STEP  SIZE' 

READ!*,*),  XSTEP 

*  — set  option  for  y  starting  t  ending  value — 
WRITE!*,*), 'DECIDE  THE  STARTING  AND  ENDING  VALIKS  OF  Y’ 
WRITE!*,*), 'INPUT  the  STARTING  VALUE  OF  Y' 

READ!*,*),  YSTART 

WRITE!*,*), 'INPUT  the  ENDING  VALUE  OF  Y' 

READ!*,*),  YEND 

*  — set  option  for  nunber  of  points  to  plot — 
WRITE!*,*), 'DECIDE  THE  NUMBER  OF  POINTS  TO  PLOT' 
WRITE!*,*), 'INPUT  THE  NUMBER  OF  POINTS' 

READ!*,*),  NPNTS 

*  — set  option  for  line  style  !IMARK) — 

WRITE!*,*), 'DECIDE  THE  IMARK  STYLE  DESIRED' 

WRITE!*,*), 'INPUT  THE  IMARK  STYLE' 

READ!*,*),  IMARK 

*  ---set  option  for  grid  style--- 
WRITE!*,*), 'DECIDE  THE  GRID  STYLE  DESIRED* 

WRITE!*,*), 'INPUT  THE  X  GRID  STYLE' 

READ!*,*),  IX 

WRITE!*,*), 'INPUT  THE  Y  GRID  STYLE' 

READ!*,*),  lY 

*  — set  option  for  output  device — 

WRITE!*,*), 'YOU  WILL  NOW  DECIDE  WHERE  TO  SEND  THE  OUTPUT' 
WRITE!*,*), 'ENTER  0  FOR  OUTPUT  TO  THE  SCREEN' 

WRITE!*,*), 'ENTER  1  FOR  OUTPUT  TO  THE  LASER' 

READ!*,*),  IPROPT 
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C  CALL  P0EV('LN03',IEER) 

IF  (IPROPT.EQ.O)  THEM 
CALL  PGPX 
ELSE 

CALL  LN03I 
ENOIF 


•  ••FUNCTION  EXECUTION  TO  DETERMINE  PLOTTING  POINTS^^ 

OPEN  (20,FILE-'SI N4002SM2 . DAT • , STATUS* 'OLD • ) 

CALL  F(XARAY,YARAY,NPNTS,XSTART,XENO) 


•  •••---LEVEL  ONE  COMMANDS-— ••• 

CALL  HUROT( 'COMIC') 

CALL  PAGEdl. 0,8.5) 

CALL  NOeRDR 

CALL  AREA20(9.0,6.0) 

* 

•  •••---LEVEL  TWO  COMMANDS- --••• 

• 

CALL  SUISSB 

CALL  HEAOINCAOO  HZ  SINE  WAVES' ,100, 1 .5,2) 

CALL  HEAOIN( 'SIGNAL  LENGTH  *  0.128  SEC  /AMP  *  2.0S‘,100,1.5,2) 
CALL  SWISSM 

CALL  XNAMECTIME  (sec)S‘,100) 

CALL  YHAMEC AMPLITUDES ',100) 

CALL  GRAF(XSTART,XSTEP,XEND,YSTART, 'SCALE', TEND) 

* 

•  *•*... level  three  COMMANDS- --••* 

CALL  CURVE(XARAY.YARAY,NPNTS.IMARK) 

CALL  FRAME 
CALL  GRID(IX,IY) 

•  ***.. -closing  COMMANDS- --••• 

CALL  ENDPL(O) 

CALL  DONEPL 
END 


*  ---FUNCTION  SUBROUTINE  TO  DETERMINE  PLOTTING  POINTS- -- 

SUBROUT  I NE  F ( XARAY , YARAY , NPNTS , XST ART , XENO ) 

• 

REAL  XARAY (NPNTS ) , YARAY ( NPNTS ), XST ART , XENO , STEP, XTEMP 
INTEGER  NPNTS 

•  ---READ  IN  THE  XARAY  AND  YARAY  VALUES- -- 

DO  100  I  *  1, NPNTS 

READ<20,^)  XARAY(I),  YARAY(I) 

100  CONTINUE 
END 
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2.  3-D  Plottiag  Routioe 


PROGRAM  UDFPLTVER3 


* 


THIS  PROGRAM  USES  THE  GRAPHIC  PACKAGE  DISSPLA  TO 
PLOT  A  3-D  FUNCTION. 


***-- -DECLARING  VARIABLES---*** 

REAL  RUDF(32768) 

INTEGER  IPROPT.IXPTS.IYPTS 


*  ***---INITIALI2ING  VARIABLES  AND  SETTING  GRAPH  OPTIONS---*** 
OPEN ( IS, F I LE> ■ RSUDF .OUT STATUS- 'OLD • ) 


DO  100  I  >  1,32768 
READdS,*)  RUDF(I) 
100  CONTINUE 


call  pdev('ln03',  ieer) 

*  **FUNCTION  EXECUTION  TO  DETERMINE  PLOT** 

*  ***---LEVEL  ONE  COMMANDS---*** 

CALL  PAGE(8.5,11.0) 

CALL  AREA20(8.0,7.0) 

*  ***---LEVEL  TWO  COMMANDS---*** 

CALL  V0LM30(6.,6.,A.) 

CALL  HEADIN('200  t  600  Hz  SINE  WAVES', 100, -1.5, 2) 

CALL  HEAD I N<' FREQUENCY  RESOLUTION  ;  5  X  35$' , 100,-1 .5,2) 
CALL  X3NAME('FREOUENCYS',100) 

CALL  YSNAMEC  ',1) 

CALL  23NAME(' AMPLITUDES ',100) 

CALL  VUANGL(-90.,0.,30.) 

CALL  GRAF30(0,, 1250. /5., 1250.0, 

:  0.,. 2048/2., 0.2048, 

:  0., 'SCALE ',.01) 

* 

*  ***---LEVEL  THREE  COMMANDS---*** 

* 

CALL  SURMATLRWDF, 1,256,1, 128,0) 

CALL  FRAME 

*  ***-- -CLOSING  COMMANDS---*** 

* 

CALL  ENOPL(O) 

CALL  DONEPL 

* 

END 
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